package scu.maqiang.cte;

import java.io.FileWriter;
import java.io.IOException;
import java.util.Arrays;

import scu.maqiang.fes.BVPType;
import scu.maqiang.fes.FES3H201;
import scu.maqiang.fes.FES3H202;
import scu.maqiang.fes.FES3H20CTE;
import scu.maqiang.mesh.Mesh3H20;
import scu.maqiang.mesh.Mesh3H8;
import scu.maqiang.mesh.Tecplot;
import scu.maqiang.numeric.Direct;
import scu.maqiang.numeric.MVO;
import scu.maqiang.numeric.NewIterSSolver;
import scu.maqiang.numeric.SRMatrix;

public class BrickResultH20_DC1 {

	public static void main(String[] args) {
		double kxx = 1.0, density = 1.0, specificHeat = 1.0;
		double E = 2.0, Nu = 0.3, alpha = 0.05;
		double beta = E / (1 - 2 * Nu) * alpha;

		double theta = 0.5, T = 1.0e1, dt = T / 1000;

	    Mesh3H8 mesh3D0 = new Mesh3H8().cube(2, 2, 4);
	    mesh3D0.scale(0.5, 0.5, 1.0);


	    Mesh3H20 mesh3D = new Mesh3H20();
		mesh3D.fromMesh3DH8(mesh3D0);

		FES3H201 ffs = new FES3H201(mesh3D);
		FES3H202 fsV = new FES3H202(mesh3D);

		FES3H20CTE fs = new FES3H20CTE(mesh3D, ffs, fsV);

		SRMatrix L = new SRMatrix(fs.GetNdof());
		SRMatrix R = new SRMatrix(fs.GetNdof());

		double[] x = new double[fs.GetNdof()];
		double[] RHS = x.clone();

		double[] temp = new double[ffs.GetNdof()];
		double[][] dispUV = new double[mesh3D.getnDim()][ffs.GetNdof()];

		double[] ex = temp.clone();
		double[] ey = temp.clone();
		double[] exy = temp.clone();
		double[] ez = temp.clone();
		double[] tx = temp.clone();
		double[] ty = temp.clone();
		double[] txy = temp.clone();

		double MM = E * (1 - Nu) / ((1 + Nu) * (1 - 2 * Nu));
		double TRef = 8.0 * 35.0 / 13.0;
		double AA = MM + TRef * beta * beta / (density * specificHeat);

		double[] coef = { density, specificHeat, kxx, TRef, theta, E, Nu, beta, dt };
		fs.assembleDierctCouplingMatrix_Type1(coef, BVPType.COMMON, BVPType.COMMON, BVPType.CONSIST_MASS, BVPType.COMMON, L, R);

//	    L.sort();
//	    System.out.println(L);
//	    System.exit(0);

		fs.applyBCHeat_MBN(L, 2);
		fs.applyBCDisp_MBN(L, Direct.X, 3, 4, 5, 6);
		fs.applyBCDisp_MBN(L, Direct.Y, 3, 4, 5, 6);
		fs.applyBCDisp_MBN(L, Direct.Z, 1);

		// System.out.println(L);
		NewIterSSolver solver = new NewIterSSolver(L);

		double time = 0.0;

		int lineNum = 100;
		double dis = 1.0 / lineNum;
		double[] lineX = new double[lineNum + 1];
		double[] lineY = lineX.clone();

		Arrays.setAll(lineX, i -> i * dis);
		String dir = "BrickResultH20\\";
		String XYLineFile = dir + "Line150.dat";
		Tecplot.lineXYTime(XYLineFile, lineX, lineY, 0.0);

		String heatFile = dir + "Heat150.dat";
		mesh3D.toTecplot(heatFile, 0.0, x);

		String deformFile = dir + "Deform150.dat";
		mesh3D.toTecplotDeform(deformFile, 0.0, dispUV, 0.0);

		String expansionFile = dir + "expansion150.dat";
		String temp00File = dir + "temp00_150.dat";
		FileWriter fw1 = null;
		FileWriter fw2 = null;
		try {
			fw1 = new FileWriter(expansionFile);
			fw2 = new FileWriter(temp00File);
			fw1.write("0.0\t0.0\n");
			fw2.write("0.0\t0.0\n");
		} catch (IOException e) {
			// TODO Auto-generated catch block
			e.printStackTrace();
		}
		while (time < T) {
			time += dt;
			System.out.println("time: " + time);
			fs.applyBCHeat_MBN(RHS, 1.0, 2);
			fs.applyBCDisp_MBN(RHS, Direct.X, 0.0, 3, 4, 5, 6);
			fs.applyBCDisp_MBN(RHS, Direct.Y, 0.0, 3, 4, 5, 6);
			fs.applyBCDisp_MBN(RHS, Direct.Z, 0.0, 1);

			solver.PCGSSSOR(RHS, x, 1.5, 1);
			fs.extractUVH(x, temp, dispUV);
			System.out.println("T min: " + MVO.min(temp) + "\t max: " + MVO.max(temp));
			System.out.println("U min: " + MVO.min(dispUV[0]) + "\t max: " + MVO.max(dispUV[0]));
			System.out.println("V min: " + MVO.min(dispUV[1]) + "\t max: " + MVO.max(dispUV[1]));
			System.out.println("V min: " + MVO.min(dispUV[2]) + "\t max: " + MVO.max(dispUV[2]));
			R.mul(x, RHS);

			mesh3D.toTecplot(heatFile, time, temp);
			mesh3D.toTecplotDeform(deformFile, time, dispUV, (1 - Nu) / (alpha * (1 + Nu)));
			// Arrays.setAll(lineY, i -> ffs.getValue(temp, new double[] {0.25, 0.25,
			// lineX[i]}));
			// Tecplot.lineXYTime(XYLineFile, lineX, lineY, time);
			try {
				fw1.write(time + "\t "
						+ ffs.getValue(dispUV[2], new double[] { 0.25, 0.25, 1.0 }) * (1 - Nu) / (alpha * (1 + Nu))
						+ "\n");
				fw2.write(time + "\t " + ffs.getValue(temp, new double[] { 0.0, 0.0, 0.0 }) + "\n");
			} catch (IOException e) {
				// TODO Auto-generated catch block
				e.printStackTrace();
			}

		}
		try {
			fw1.flush();
			fw1.close();
			fw2.flush();
			fw2.close();
		} catch (IOException e) {
			e.printStackTrace();
		}
	}

}
